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Abstract — We present a novel technique by which highly-segmented 
electrostatic configurations can be solved. The Robin Hood method 
is a matrix-inversion algorithm optimized for solving high density 
boundary element method (BEM) problems. We illustrate the 
capabilities of this solver by studying two distinct geometry scales: 
(a) the electrostatic potential of a large volume beta-detector and (b) 
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the field enhancement present at surface of electrode nano-structures. 
Geometries with elements numbering in the O(IO^) are easily modeled 
and solved without loss of accuracy. The technique has recently been 
expanded so as to include dielectrics and magnetic materials. 

1. INTRODUCTION 

Numerical modeling of the electromagnetic properties of complex 
systems is often a necessity for many scientific and engineering 
disciplines. Knowledge of the electric and magnetic fields often 
provides valuable insight into the system's performance, optimization 
parameters, and inherent response. As these systems grow in size 
and complexity, so do the demands upon the model employed. Even 
within a well-understood framework such as Maxwell's electromagnetic 
theory, the user is often confronted with a choice between accurate, 
model-independent solutions or timely results. 

The development of computational power in the last decades has 
resulted in a large number of available efficient methods for the solution 
of these potential problems in electrostatics, as well as in other areas 
of physical, engineering and technological applications. According to 
the basic approach towards solving electrostatic problems, methods 
can be roughly divided into either finite element methods (FEM), 
boundary element methods (BEM) and finite difference (FD) methods. 
BEM methods in particular have risen in wide-spread use over recent 
decades [D El 13 |ll|5l [6] , as it often provides the most accurate means to 
describe a given system. Employing BEM from a practical standpoint, 
however, becomes somewhat difficult for highly segmented geometries, 
as the requirements on computational memory grows quadratically 
with the number of sub-elements. 

In this article, we describe how a new BEM-solving algorithm, 
called Robin Hood, elegantly solves the memory problem, providing 
accurate results with reasonable computation times. The technique 
was first introduced in 2006 ff, "8] , and was recently expanded to include 
dielectric and magnetic materials. We illustrate its applicability and 
scope by studying both micro- and macro- scale systems. Section [2] 
describes the nature of the problem to be solved. Section [3] describes 
the Robin Hood method and its basic performance characteristics. 
Section |4] describes its use in solving real- world systems. We conclude 
in Section [5] with a summary and outlook for future applications. 
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2. THE ELECTROSTATIC BOUNDARY VALUE 
PROBLEM 

The Boundary Element Method (BEM) is a computational technique 
for solving linear partial differential equations. Compared to other 
popular methods (such as the Finite Element and Finite Difference 
Methods j^) designed to accomplish the same goal, the BEM differs 
in many respects, favoring its use in an important subset of problems. 
Instead of discretizing the entire region of interest, the main technique 
of the BEM is to discretize only the surfaces of the geometries in the 
region. This effectively reduces the dimensionality of the problem and 
facilitates the calculation of fields for regions that extend out to infinity 
(rather than restricting computation to a finite region) [10]. These 
two features make the BEM faster and more versatile than competing 
methods when it is applicable. 




Figure 1. A two-dimensional graphical depiction of the regions O (in 
white) and S (in grey) (see Fig. [T]). The vector n describes the unit 
normal to the boundaries (F^ and Fs) of these regions, and x is the 
observation point. 

We begin by defining a two-dimensional inner region fi, bounded 
by a piecewise smooth contour F^ with clockwise orientation. The 
region $7 is bounded by an external region E, which is bounded by 
a circular, counter-clockwise contour Fs of radius i?s and internally 
by the contour F^. To derive the principle equation of the BEM, one 
begins by applying Green's second identity to a region S encapsulated 
by mixed (Dirichlet and Neumann) boundaries Tq and Fs, for two 
arbitrary scalar functions U and W 
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J (UV'^W-WV^U^ dA' = 

•^E+rn dn ^ dn) ' 

We define the electric potential, U (x) in Eq. [T] such that the 
Laplace relation is satisfied, i.e.: 

V^U{x) = 0. (2) 

Defining W as the fundamental solution to the Laplace equation, 
we find: 

W = G{x,x') = -^—^. (3) 
Aireo \x - x'\ 

When extending Fs out to infinity and performing contour 
integrations about the boundaries of S, one is left with the relation: 



c,ix).Uix)= / U—-G—)dS', (4) 




ci{x) = { i-'^ xGTn , (5) 



where 6q is a boundary angle of integration whose pivot is located at 
Fa (see FigM. 

Often, Eq. |4] is numerically solved by discretizing the boundary 
Fq (this technique is known as the direct BEM). Another technique, 
known as the indirect BEA^^ applies the above formalism to the inner 
region Q (as in Fig. [s]), yielding a similar relation to Eq. [4| 



C2(x) • U{x) = I U^-G—] dS', (6) 




C2{x) = { t x€Fn , (7) 



which differs from Equation [4] by the definition of the boundary angle, 
the boundary's orientation and the direction of the boundary normal. 



JThis technique is also known as the Source Element Method (SEM), Source 
Integration Method (SIM), and Charge-Density (Integral) Method. 
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Figure 2. Same as Fig. [T] except the F^ boundary is deformed to 
include a small circular segment of radius e to avoid the singularity of 
V^M^(r) at the boundary. In the limit e — t- 0, the original boundary is 
recovered. 




Figure 3. The equivalent figure to Fig. [T] for the internal region Q. 
The shaded region denotes the domain of interest for computation. 

By reversing the boundary normal in Eq. [6] and restricting the domain 
to S, Eq. [6] simplifies to 

By matching the boundary conditions and summing Eqs. |4] and 
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|8j one finds 




(9) 



Defining a to be the sum of tlie fluxes across Tq: 




(10) 



U{x) = f G{x,x')a{x')dS' 



(11) 




Figure 4. Discretization of Tq into n sub-elements for numerical 
computation. Charge densities are constant along a sub-element. 

If the distribution of charges is known apriori, such a calculation 
would be relatively trivial to solve. Of course, the distribution of 
charges is not known, but needs to be solved for the given geometry. 
What is typically known, however, are the boundary conditions of 
the system. In the case of conducting surfaces, this represents the 
electrostatic potential of the surface of the given conductor. Even 
when the actual value of the potential is not necessarily known, such 
as in the case of the isolated conductor, the condition of equipotential 
surfaces still holds (i.e. the potential is constant across the entire 
surface). As will be shown later, the case is slightly altered when 
dealing with dielectric materials, though the mathematical foundations 
are unchanged. 

In solving for the distribution of charges, the potential is therefore 
evaluated at the surface boundaries. One often resorts to numerical 
methods involving discretization of the surface in order to solve for 
the charge distribution and hence determine the potential at any point 
in space. The discretization is often carried to the point where it 
is possible to treat the surface charge density for a given element as 
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uniform over the sub-element surface dS' . In this case, Eq. 11 can be 
re- written in matrix-like notation: 

N 

Ui = Y,hj-aj, (12) 
i=i 

where Ui is the potential of the sub-element i, aj is the charge density 
of sub-element j of total elements, and lij is given by the expression: 

h = J^I T^r (13) 
47reo J AS J \ri - rj] 

The element lij therefore represents Green's function correlating 
sub-element i at to sub-element j at rj . Inversion of Eq. [T2] yie lds the 
relevant discrete charge distributions from which Equation |11| can be 
fully evaluated. However, for geometries where extremely large number 
of elements exists (warranted either by the required accuracy, the 



extent of the surface involved, or both), numerical inversion of Eq. 12 
can impose often severe requirements, particularly on convergence ana 
memory. As we shall see in Section |3.1[ the Robin Hood method 
allows for elegant solution to the inversion/memory problem. Before 
we leave the section, we expand the formalism so as to also include 
linear dielectric and magnetic materials 



2.1. Extension to Dielectrics and Magnetic Materials 

In considering the calculation of electrostatic potentials in the presence 
of dielectric media, one no longer can impose the condition that each 
surface can be treated as an equipotential. Nevertheless, it is certainly 
possible to recast the problem in terms of a matrix equation where one 
solves for the induced charges along the boundaries of the system [12j. 
We still consider the case that Laplace's equation on the potential still 
holds: 



V^U{x) = 0. (14) 

Let us take the boundary condition that at the insulator-insulator 
interface the displacement vector D is continuous across the surface. If 
the surface represents the boundary between two dielectric materials, 
then the boundary condition on the surface imposed by Maxwell's 
equations results in: 

etEt■h^-er-Er.h^ = 0, (15) 
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where ef is the permittivity above and below the surface of the 
subclcmcnt i, is the value of the electric field at sub-clcmcnt i, and 
■hi is the surface normal vector at the insulator-insulator interface. The 
boundary condition can be re-expressed also in terms of the gradient 
of the potential: 



[rii 



V)C/+-e-(n,.V)t/-=0. (16) 



Here, the expression [hi ■ V) represents the gradient oriented 
along the direction normal to the surface i. In order to evaluate the 
electric field at a selected point, it is convenient once again to express 
everything in terms of their integral forms. In this case, we can take 
advantage of the kernel evaluation of Green's function: 

(n, . V)G(r- ,f,) = (17) 
47reo — '^jl 

Hence, the electric field can be expressed as the following sum over 
finite elements: 



With some manipulation of the above equation, we can re-express 
it in matrix notation, just as we did for the scalar potential: 

N 

where for i ^ j we have 



Vij = / |!' !^|3 • fiidSj, (20) 

47reo JASj In -rjf 

and for i = j we have 



1 Mi±iO pi) 



47ren 



In the case of dielectrics, the boundary condition forces the 
constraint vector, ^i, to be zero for all elements. The matrix equation 
can easily be extended to include a mixture of dielectric and conducting 
surfaces, as long as one calculates all contributions to the electric 
field and electric potential when evaluating the appropriate matrix 
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elements. Formulated in this manner, and taking advantage of the 
linearity of electric fields, one can impose the same algorithm to solve 
conductors and insulators. 

The applicability of the Robin Hood method to both Neumann 
and Dirichlet boundary conditions implies that one can extend the 
method to the calculation of magnetic fields, even in the presence 
of materials with magnetic permeability [13]. For the case of 
magnetostatics in conductive media (no magnetic materials), one 

wishes to compute the magnetic field B{x) due to specified current 
configuration. It is best to describe the system in terms of the magnetic 
vector potential. A, where B = V x A. If we decide to work in the 
Coulomb gauge 



where K is the surface current density and is the permeability of 
free space. If the current density is zero within the volume considered, 
then it is possible to express the magnetic field in terms of a magnetic 
scalar potential. In such cases, one often has the knowledge of the 
current being distributed to the system, hence, performing a matrix 
inversion to solve for the free currents is not necessary. However, once 
magnetic materials are introduced into the system, one no longer has 
knowledge of the induced currents, and one therefore returns to the 
same problem as posed in the case of dielectric materials. Fortunately, 
we use an approach similar to that employed for dielectric materials to 
solve for such cases. 

Let us define the magnetic field above and below the boundary 
of some magnetic material as B^. Within the Coulomb gauge, the 
magnetic field can be expressed in terms of the vector potential A: 




then the vector potential can be written in a very similar manner as 
in electrostatics: 




(23) 



N 




(24) 
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n, x(^B,+ -— Br) = 0. (25) 
Substituting Eq. [24] into Eq. [25] yields the condition: 

X (ip - 4) E(V X A,(rl)) + (i^ + = 0. (26) 

Dividing out the common term, one finds: 

Y: n.x(VxA,(r-;)) + y ^^;/ k. = 0. (27) 

The magnetic vector potential term can be simplified further by 
realizing that the vector A is continuous along the boundary 

nxB = nx(VxA) = V(n • A) - (n • V)A. (28) 

In order to simplify the constraints on each of the components of 
the current K, let us consider both the vector normal to the surface 
o , Tii 5 as well as the two independent components tangential to the 
surface, i\, where I = 1,2. The boundary conditions impose different 

conditions on the tangential and normal components of K. Specifically, 
the tangential boundary conditions can now be written as: 



E ^ / '^^^^TTli • (*'(^^ • - n.(tl • K,)) (29) 



j=l,jf^i 



+ ^,f-'^^/,^ ^l-Ki) = 0. (30) 

If once again we assume the surface current does not change across 
the integration, we can express this as a matrix whose elements are 
given as follows 

N 

= E X^J ■ {nS\ ■ K,) - i\{h, ■ Kj)), (31) 
where for i ^ j we have 

Ho f (ri-rj) 
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and for i = j we have 

& = (f5)?^^<f±ii>n. (33) 

Since at the surface i the current is constrained to flow along the 
surface, there is no normal component • Kj. 

These conditions completely define induced surface currents and 
surface charges from all elements. For linear materials, these equations 
can be solved to obtain the values of a and K across all surfaces by 
means of once again inverting an x matrix. Note that we have 
restricted ourselves to the case where there are no external source 
charges or currents in the volume. Had this assumption been relaxed, 
the boundary vector *j would take on values other than zero. 



3. TECHNIQUE AND IMPLEMENTATION 
3.1. The Robin Hood Method 

Let us consider the case where an isolated conducting sphere comes 
in close proximity to a charge (or collection of charges) some distance 
away. Although the total charge on the sphere is zero, the charges 
rearrange themselves so as to ensure that the electric potential 
everywhere on the conduct is constant, thereby ensuring Gauss' law 
is satisfied. If there is an imbalance in the potential, charges move 
and swap until the proper potential at the surface is reached. Nature, 
of course, accomplishes this nearly instantaneously but the rules that 
guide the distribution of charges can be simplified such that computer 
algorithms can be made to mimic the effect. 

The Robin Hood method optimizes this natural strategy to 
conform on how computer algorithms are designed, while maintaining 
the simplicity of the rules illustrated briefly above. We will discuss the 
details of the Robin Hood implementation in the next section. 



3.2. Implementation 

Consider once again our example of the isolated conducting sphere in 
close proximity of a single charge of value +q located at some distance 
away. As flrst step, which embeds the approximation scheme used 
in this approach, we subdivide the geometry of our large sphere into 
A'' discrete triangles and assume a constant charge density across the 
surface of each triangle. The discretization determines the level of 
accuracy achievable by the model. Since the charge distribution is 
unknown, we randomize the individual charge distributions for each 
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triangle. If the object had some finite charge, the constraint would 
be imposed that the sum of all charge add up to the assigned charge. 
The potential is then calculated from this initial assignment of charge 



distributions, in accordance with Equation 12 



The bulk of the computing time is spent in calculating the matrix 
coefficients lij, which essentially reduces to computing the potential 
from a triangle at some arbitrary point Xi. We have the option of 
employing a variety of techniques in computing this potential. In cases 
where accuracy precedes timing needs, we have the option of computing 
the potential from a triangle analytically. Such a computation is a bit 
more time consuming, as the evaluation has terms that depend on 
sinh~^(x) and tanh~ (x) functions. Alternatively, one can use self- 
refining mesh and simply compute the multipole expansion up to the 
quadrupole term to obtain the value of the potential at the given point. 
The criterion for accuracy using this approach is expressed as a ratio 
(let us denote it by v) of the distance \xi—Xj \ and the typical size of the 
right-angled triangle which is characterized by the length of the larger 
leg of the triangle. When the ratio v is larger than some pre-determined 
cutoff (in our case, 27), the monopole term alone is enough to obtain 
the value of the potential at point Xj with 4 digit accuracy. If the ratio 
v is smaller that 27 but larger than 5.7 then monopole and quadrupole 
term together yield a 4 digit accuracy. If the ratio u is smaller than 5.7 
then the initial triangle is to be divided into four new triangles. In that 
case, the potential becomes the sum of 4 contributions from the new 
triangles, which are considered homogeneously charged bearing the 
same charge density as the initial triangle. A full list of the expressions 
used in the triangle potential evaluation can be found in the Appendix. 

Once the potential is evaluated at each boundary element in the 
problem, the algorithm then performs a search for the two elements 
that differ most from the equipotential surface. Assume that two of the 
worst triangles are electrodes m and n. The equipotential condition 
imposes that these two elements have the same potential. Therefore, a 
small amount of charge is moved from one electrode to the other such 
that the equipotential condition is met; that is, U!^ = U^. The new 
potential due to this redistribution is therefore given by the expression 



U!^ = Um + Imm ' + Imn ' ScTn, (34) 

Un = Un + Inn ' ^(^n + Inm ' ^(^ni- (35) 

where Um,n is the potential at the two selected points. For the case of 
two sub-elements, there is an exact solution to the amount of charge 
that needs to be moved from one element to the other. 
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e An{Un Urn) fna\ 

^''"^ = TTi — -I \ + A (I -I y ^^^^ 

" " ~ATi -/ \ + A (I -I V ^ ' 

where 5am,n is the charge density exchanged and Am,n is the area of 
the triangle sub-element. Note that the sum of the two charges is zero; 
i.e. Sqtotai = df^m • A„i + 5(Jn ■ An = 0. Though we have used the 
example of an isolated sphere to illustrate the essential elements of the 
algorithm, we can certainly extend the case for a conductor held at a 
fixed potential. When the sphere is held at a fixed potential, Uq, the 
shift of charges between the two electrodes is altered to the following 
expression: 



{Uq — Um)Inn — {Uq — Un)Imn ,oo\ 

(>(^m = f f — ^7 — f , (38) 

^ {Uq — Un)Imm — {Uq — Um)Inm yr,n\ 
S(T„ = ^ _ — ^ . (39) 

Having solved for the shift in charge, the potential for every 
sub-element is subsequently updated to reflect the new charge 
configuration. The method is said to converge once the maximum and 
minimum values of the potential on an individual electrode is below 
some user-defined value. In all examples used to date, convergence 
appears to be scale- independent, i.e. it continues to fall exponentially 
until machine precision is achieved. As can be seen from the above 
description, the algorithm is extremely simple in implementation, yet 
encodes the basic natural properties imposed by Gauss' law. 



3.3. Some Generalizations 

So far we have restricted ourselves to the case where just two electrodes 
are involved in the charge exchange. We can generalize the method to 
the case of M electrodes exchanging charge, where 1 < M < A?". For 
the case where the electrodes are held at fixed potential, the charge 
exchange solution is equivalent to the inversion of a M x M matrix: 

M 

Ul = Ui + Y,IijSaj, (40) 

where U- is the adjusted potential after a small amount of charge 
Sqj = ASj ■ Saj is exchanged between M selected electrodes. Both 
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indexes i and j range from 1 to M. Note that for the case that 
M = 2 and the target potential is known (i.e. U- = Uq for all selected 



electrodes), one reproduces Eq. 39 



The same generalization can be made for the insulated charge case. 
Here, one must explicitly introduce the conservation of charge as part 
of the matrix to be solved. In general, one can impose the conservation 
of charge as: 



M 



J25q,=0. (41) 



As for the condition imposed by the equipotential surface, all the 
potentials at each altered electrode are identical. This condition can 
be enforced via the recursive relation: 

M 

- Ii,j)6a, = U,- Ui+i (42) 

The index i ranges from 1 to M — 1. The combination of the 
two constraints produces an M x M matrix equation which can be 
inverted to solve for Uj. Note that for the isolated conductor we 
require that M > 1. In the case where M = 2, one reproduces 



Eq. [36) It is interesting to study the properties of these solutions 
for various limits on the value of M. For the case M — )• 1, the Robin 
Hood method becomes a variant of the Gauss-Seidel method, albeit 
it is implemented with a much lighter memory footprint. In the limit 
M — )• A^, it approaches the original matrix inversion problem. In 
general, the convergence behavior is best for M — )■ 1. 

There is one last generalization that can be useful in practical 
applications. In most study cases, one is not just interested in a 
single configuration, but in the behavior of the system under different 
boundary conditions. As long as the geometry remains fixed, it 
is possible to solve the charge configuration to allow any arbitrary 
combination of voltages. 

Consider the case where one can break up the boundary into 
different groups of electrodes (each group corresponds to a subspace of 
Ffj), each of which is held at a distinct potential that is common to all 
the elements in that subgroup. It is now possible to solve for the charge 
configuration for the case where all the elements in a given subgroup 
are held at a common potential, ^q, while all the other potentials are 
held at ground; that is: 
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Ut = J2^^^^h (43) 

j 

where k denotes the specific group, Ui is zero everywhere except for 
the electrodes in group k, and represents the charge solution for 
the group. This procedure is repeated for each and every group in the 
set. Because the electric potential everywhere in space depends on the 
linear superposition of the charges, we can now construct the electric 
potential due to any arbitrary setting of the potentials, on each of 
these groups: 

= E (44) 

k=l 

where Ng is the total number of electrode groups in the configuration 
and a, is a scalar defined as 



ai = ^i/^o. (45) 

Although the computation time obviously increases approximately 
by Ng, further computation is no longer necessary even if the potentials 
on the electrodes are altered. 



3.4. Convergence Properties 

The Robin Hood approach is expected to converge as long as the 
following three conditions are satisfied: 

(i) The system is linear. 

(ii) The matrix is diagonally dominant. 

(iii) The charge configuration corresponding to the boundary condi- 
tions is unique. 

Conditions (i) and (ii) are satisfied by the properties of Maxwell's 
equations and Green's function. Indeed, since the potential falls like 
1^ j*^^ I , the self-term diagonal element in the matrix is always greater 

than its neighbors, and so convergence is guaranteed. One could argue 



that in expression [36] the denominator Imm + Inn — Imn — Inm could 
become zero thus causing the Robin Hood approach to fail. However, 
for a homogeneously charged triangle the potential has the largest 
value at the barycenter of the triangle, which implies that for any 
other element n the relation > I^m always holds. In particular, if 
it happens that /„„ is equal to Inm, it implies that the two different 
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triangles share the same barycenter, which means that they are either 
intersecting or completely overlapping. Neither of these situations 
should be allowed because in that case the whole electrostatic problem 
becomes under defined thus allowing for infinitely many solutions. 
Therefore, in well defined geometries, condition (iii), and by inference 
condition (ii), are always satisfied. 

3.5. Advantages 

Despite its extremely simple algorithm, the Robin Hood method offers 
a number of concrete advantages over other matrix inversion methods. 
The main advantage of the technique is its memory footprint. Under 
the Robin Hood approach, there is no need to store individual matrix 
elements during the computation. Rather, only the potential at each 
element is stored, while the matrix elements are simply re-calculated 
when updating the potential due to the charge exchange. As such, the 
memory requirements grow linearly with N, instead of quadratically. 
This feature allows for solving the charge configuration even for 
extremely segmented geometric configurations previously avoided. 

The other advantage inherent to the technique is that it is 
extremely parallel in nature. The computation of potential can 
be distributed to any number of CPUs, since each computation is 
essentially isolated. In fact, with the exception of the charge exchange 
calculation itself, the entire sequence can be passed to a cluster of 
machines to facilitate computation. For this paper, the calculations 
are typically performed on clusters using the Message-Passing Interface 
(MPI) protocol with great success in reducing real calculation time. 
The Robin Hood method also has recently been extended to function 
on Graphical Processing Units using the OpenCL standard pS], and 
has been verified on both ATI and NVIDIA graphics cards. 

One last advantage that Robin Hood provides is that the 
computation time grows quadratically with A^, while using direct 
solvers the growth is cubic. 

4. CASE STUDIES 

4.1. Verification Examples 

The Robin Hood algorithm has already been benchmarked in previous 
publications against standard examples used in electrostatics. Most 
notably, the Robin Hood method has been used to calculate the 
capacitance of a cube. The largest cube calculated was represented by 
202,800 triangles, obtaining a capacitance value of 0.66067786±8x 10~^ 
in units of l/ineQ. The calculation represents one of the most accurate 
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values of the capacitance of the cube calculated to date [3 [8]. The 
algorithm is not limited in terms of the number of triangles employed 
and, when used in combination with parallel processing, can be 
remarkably fast and efficient. 




Figure 5. Left: Two nested conducting spheres, one held at ground 
potential while the other held at a fixed potential of 10 V. Right: 
Distribution of potential value for (solved) nested sphere geometry on 
left. The potential distribution is sampled throughout the volume of 
the innermost sphere. 

The capacitance is a poor metric for testing the accuracy of 
the method, in part because it is an integral measure of the system 
variation. A more compelling comparison is one that tests certain 
fundamental principles, such as Gauss' law. We consider the case where 
we have a set of nested metal spheres, where the inner sphere is held at 
a fixed potential (10 V) while the outer shell is held at ground. Gauss' 
law predicts that the electric potential inside the innermost sphere is 
fixed, regardless of whether the inner sphere is concentric with the 
outer sphere or not. For our study, the sphere was subdivided into 
20,000 triangles using GMSH's Delaunay algorithm [38j . The accuracy 
of the calculation is determined by the metric 

N 

e = ^|C/(x,)-f/target|. (46) 

i 

The quantity e is sampled throughout the volume inside the 
innermost sphere. Results for our test can be seen in Figure [5j For the 
example cited above, the average voltage within the sphere was within 
±15//V of its expectation value. 
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4.2. M- and N- Dependence Studies 

In order to test some of the fundamental properties of the Robin 
Hood method, we made electric field computations of a simple dipole 
electrode system. This consists of two electrodes having opposite 
potentials (1000 V and -1000 V), each of them having a half-cylindrical 
shape, with some gap between them (see Figure [6]). We discretized the 
surface of these electrodes with independent triangle surfaces. All 
iterations were started with zero charge density on all triangles. 




-l,857e-6 1.857e-6 



Figure 6. A model of the dipole configuration used for M- and N- 
dependence studies. The color indicates the surface charge density 
distribution of the object. 

We investigated the M-dependence of the potential accuracy and 
of the iteration number, with a discretization of = 3600. Table [T] 
contains results of the dipole charge solution computed with 6 different 
M values. We study both the relative potential accuracy Aj-^i values 
achieved by the Robin Hood method with fixed M ■ n product and with 
fixed relative accuracy. The relative accuracy Arei denotes the maximal 
deviation between the input potential and the simulated potential, over 
all triangle subelements. Since the computation time t of the Robin 
Hood method is proportional to M-n, the computation time is roughly 
constant. As Table [l] shows, we found the best accuracy with M = 1 
(the Gauss-Seidel limit of the Robin Hood method). For a fixed relative 
accuracy, it appears that M = 1 is also the best choice for minimizing 
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the computation time. 

Table [2] presents the number of iterations n for various accuracy 
levels, using the same discretization as above, and M = 1. It seems 
that n is approximately a linear function of log^rei- Figure [7] shows 
the time dependence of the RH solver as a function of the number of 
sub-elements once the relative accuracy is fixed. This time evolution is 
nearly quadratic in N. These results are also shown in Table [3} As a 
testament to the vast improvement to be gained from GPU processing, 
Figure [7| also shows the time evolution for the dipole configuration as 
solved with GPU implementation. 

We have repeated these simulations with several other electrode 
systems: unit square and unit cube, various electrodes containing both 
rectangles and wires as subelements etc. In all cases, we obtained the 
same results as above: 

• The Robin Hood method is fastest with M = 1 (Gauss-Seidel 
limit). 

• The iteration number increases logarithmically with the relative 
potential accuracy. 

• The iteration number increases linearly and the computation time 
quadratically with the number of the discretization subelements 
(for the CPU implementation). 



M 


1 


2 


3 


4 


6 


10 


^rcl 


7.23 • 10"^^ 


4.5 • 10"^^ 


8.3 • 10"^^ 


1.6 • 


3.9 • 


3.2 • 10-y 


M ■ n 


20999 


22580 


23547 


24532 


26796 


29750 



Table 1. Robin Hood results with a dipole electrode (A^ = 3600), 
for 6 various M values. Second row: relative accuracy A^^x for fixed 
M ■ n = 30000 value {n: number of iterations). Third row: M ■ n for 
fixed relative accuracy A^^i = 10~^. 



^rel 


10"^ 


10-4 


10"'^ 


io-« 




n 


3425 


9108 


15047 20999 


26973 



Table 2. Robin Hood results with a dipole electrode {N = 3600): 
number of iterations n for various relative accuracy values ^rei- 
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Convergence time: CPU vs GPU 




2000 4000 6000 8000 100001200014000 
# of sub-elements 



Figure 7. A plot of the time dependence versus number of electrodes 
for the dipole configuration. The plot shows the computation time 
for CPU processing (blue) and GPU processing (green) versus number 
of sub-elements. In all cases, the required relative accuracy of the 
potential at the surface of each triangle electrode is set to be 10~*. 



N 


900 


1444 


3600 


7056 


13456 


n 


5279 


8496 


20999 


40859 


77030 


n/N 


5.87 


5.88 


5.83 


5.79 


5.72 


t/N'^ 


3.0- 10-^ 


3.0 • 10"^ 


2.9-10-^ 


2.9-10-^ 


3.0 • 10"^ 



Table 3. Robin Hood results with a dipole electrode for 5 various 
discretizations, with fixed relative accuracy ^rei = 10~^, as computed 
with CPU implementation. Number of independent subelements: N , 
number of iterations: re, computation time (in seconds): t. 
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Having verified the efficacy of the method and its optimal 
performance, we now choose two very distinct examples -a large- 
scale system and a micro-scale system - to illustrate the technique's 
usefulness and versatility in real-world applications. 

4.3. Case 1: Large-Scale System: the KATRIN Focal Plane 
Detector 

The KArlsruhe TRItium Neutrino experiment (KATRIN) combines 
an ultra-luminous molecular tritium source with an integrating high- 
resolution spectrometer to gain sensitivity to the absolute mass scale 
of neutrinos. The projected sensitivity of the experiment on the 
neutrino mass scale is 200 meV at 90% C.L. |19J. A brief description of 
the overall experimental setup, currently under construction, is given 
below. 

The KATRIN experiment is based on technology developed by 
the Mainz and Troitsk tritium beta decay experiments. These 
experiments used a MAC-E-Filter (Magnetic Adiabatic Collimation 
combined with an Electrostatic Filter) [18]. The technique, which 
combines high luminosity with high energy resolution is as follows: 
Low-pressure tritium gas is stored in a high magnetic field tube. Beta- 
decay electrons from the gas, though emitted isotropically, are guided 
longitudinally along the magnetic field towards a large spectrometer, 
while neutral gas is returned to the source tube by pumps (see 
Fig. [s]). On their way into the center of the spectrometer, the electrons 
follow the magnetic field while its strength drops by many orders 
of magnitude. The magnetic gradient force transforms most of the 
cyclotron energy Ej^ into longitudinal motion. The resulting parallel 
beam of electrons runs against an electrostatic potential formed by 
cylindrical electrodes. Electrons with enough energy to pass the 
electrostatic barrier are reaccelerated and collimated onto a detector, 
while lower energy electrons are reflected. Therefore, the spectrometer 
acts as an integrating high-energy pass filter, with which we measure 
the P spectrum in an integrating mode. More detailed descriptions of 
each component can also be found in the KATRIN Design Report [T9] . 

In this study, we concentrate on the modeling of the focal plane 
detector (FPD), where the /3-decay electrons are focused and detected. 
An illustration of the FPD is shown in Figure [9} Electrons enter left of 
the first (pinch) magnet and are magnetically collimated onto a silicon 
pin diode detector. A horn-shaped electrode with a quartz insulator 
surrounds the magnetic flux lines in order to provide additional 
acceleration to electrons entering the system as a means to mitigate 
against low energy backgrounds. Calibration and pumping ports can 
be seen at the outer edges of the vacuum tube, making the system's 
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Figure 8. Schematic overview of the KATRIN experimental 
setup: (a) rear section, (b) tritium source, (c) diffrential pumping 
section, (d) cryogenic pumping section, (e) pre-spectrometer, (f) main 
spectrometer, (g) detector system. The overall setup has a length of 
about 70 m. The Robin Hood model is based on the detector system 
(g) only. 

geometry highly asymmetric. 




Figure 9. Left: the primary components of the detector system. 
Right: The Robin Hood model of the detector region, solved for the 
potentials applied to the electrode surfaces. 

The system's asymmetric geometry forces one to model the system 
using triangle discretization. The experiment is sensitive to the 
presence of low level Penning traps (instances where magnetic field lines 
stretch across two separated electrodes, both at non-zero potential). 
Such Penning traps are often a source of electrons and ions that 
contribute to the background of the experiment. Finding such traps 
requires high level of geometric discretization and accuracy. The 
memory requirement imposed by such high level discretization is often 
prohibitive for other electrostatic solvers. 

The focal plane region was broken up into several electrode groups 
(each group is defined by whether it is a dielectric or held at a specified 
potential) and sub-sequentially divided into a total of 444,821 triangle 
sub-elements of varying scale. Scaling is usually determined by one's 
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proximity to intersecting structures or to sharp edges. Calculation 
of the surface charge densities for this electrostatic configuration was 
carried out at the University of North Carolina computing grid, which 
consists of 520 blade servers, each with 2 quad-core 2.3 GHz Intel 
EM64T processors, 2x4M L2 cache (Model E5345/Clovertown), and 
12 GB of memory. The calculation was distributed using MPI protocol 
over approximately 100 CPUs. A total of ~ 1.5 x 10^ iterations were 
carried out to each a minimum relative accuracy of e < 10"*^ on all 
sub-elements. The total computation time was approximately 16 CPU- 
months (~ 4 days real time). 





CD/ICDI*log(l + ICD*l.el6l) 
20 ,0 , 20 
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Figure 10. Charge density distribution for section of detector region. 



Figure 10 shows the solved charge density distribution of 
the detector system after convergence. Despite its high level of 
discretization, presence of dielectrics, and geometric complexity, the 
charge configuration is well defined throughout the volume, and its 
convergence is achieved in a reasonable time frame. The convergence 
scale 
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as a function of iteration is shown in Figure 
charge configurations can be used to calculate various aspects of the 
electrostatic configuration, such as the presence of Penning traps and 
sources of potential vacuum discharge. The results of this study are 
the topic for an upcoming paper. 



4.4. Case 2: Micro-Scale System: Nanostructuring at 
surface 

Field enhancement on nanostructures at the surface of electrodes is of 
great practical importance, particularly in the fields of nanoplasmonics 
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x10^ 



Iteration number 

Figure 11. Relative accuracy of detector electrode groups 
(represented as different curves) as a function of iteration. The 
different curves represent different component groups which share the 
same target potential within the detector. The convergence of the 
geometry scales exponentially with the number of iterations. The 
calculation was performed until all electrode groups achieve a relative 
accuracy of 10~^. 

|24j and the development of ion emitters P2] ■ Nanoplasmonics is a bit 
more complex, since in principle it includes the time-dependence of the 
electric fields [25]. so we chose for this study to concentrate solely on 
ion emitters. 

Ion emitters are designed to produce a relatively large current by 
ionizing atoms close to the surface of an electrode. Field ionization 
occurs via electron tunneling on through a potential barrier between 
a neutral atom and an electrode maintained at a fixed voltage. The 
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probability of tunneling depends of the strength of the electric field near 
the electrode and, as such, depends on both the potential difference 
applied to the electrode as well as the specific geometry of the electrode 
itself [56]. Due to the dependence of the electric field strength on the 
geometry, modeling of the properties of the ion emitter, such as field 
enhancement, is non-trivial. In many applications, nanostructures are 
deliberately produced to maximize field-enhancement |22j . This allows 
for simpler power requirements and often enables such devices to be 
much more portable [28l |27] . 

The optimization of different geometries for increased electric 
field enhancement was studied mostly experimentally, where only a 
few vague empirical rules were derived [29] . We are aware of only 
one theoretical study of field enhancement for various geometries [2D] . 
Part of the reason for the lack of a solid theoretical framework is due 
to the absence of solvers that can provide accurate solutions of the 
electric field for relatively complex geometries with a large number 
of boundary elements. In techniques such as finite element or finite 
difference, the problem is numerically prohibitive, while for most BEM 
implementations, the memory requirement is often too large. The 
Robin Hood techniques offers a possible solution to this modeling 
problem. 

Often the field-emission electrodes used in ion emitters are covered 
with extremely thin objects, such as carbon nano-tubes. Due to 
the presence of such very sharp tips, the electric field is extremely 
enhanced compared to a smooth electrode, up to factors of 10^. Such 
enhancement can be tempered by field screening caused when other 
high field regions are brought to its vicinity, which is typically the 
case in producing such patterned electrodes. A greater density of such 
surface deformations stems from the need of large emission currents. 
Therefore, there exists a tradeoff between the field enhancement on 
the tip and number of tips per surface area [211 EO]. To study the 
relationship between field and current enhancement, we explore two 
possible geometries. 

4.5. Egg-Carton Geometry 

Field enhancement on corrugates surfaces is commonly used today in 
surface-enhanced Raman scattering (SERS) [3T], where such effects 
can increase the measured Raman signal up to 11 orders of magnitude 
due to the E'^ dependence on the electric field strength. SERS is in 
principle a time-dependent phenomenon, but for certain cases, a quasi- 
static calculation is sufficient to estimate the quality of the corrugated 
surface for the enhancement of the signal. Some other aspects of 
field enhancement on the corrugated surface can include unusual and 
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unexpected phenomena, such as X-ray emission [33l [32] . Because of 
the typical nano-structured geometry that appear in some of these 
appHcations, we have chosen a simple egg-carton structure. The surface 
is modeled by a | cos (Ar)| function where the amplitude is constant but 
the spacing -controlled by A- is allowed to vary. An example of such 



a geometry is shown in Fig. 12 




Figure 12. Depiction of the egg-carton geometry. The realistic grid of 
boundary elements (triangles) that is used in the calculation is visible. 
The peak-to-spacing ratio shown for this figure is 4 to 1. 



We consider the geometry from Fig. 12 to be an ideal metal 
electrode held at a fixed voltage. Using the Robin Hood technique 
outlined earlier, we find a solution that satisfies the boundary 
conditions at the surface of the electrode at a relative accuracy of 
e ~ 10~^. Having solved for the surface charge densities, we can then 
calculate cross-section of the potential and electric field, as is shown 



in Figures 13 and Fig. 14 Using the cross-section electric field values, 
we determine the maximum electric field strength relative to where 
the field is homogeneous. From this ratio, we can determine the field 
enhancement factor for a particular geometry. We have calculated field 
enhancement factor for different amplitude versus spacing ratios and 
results are shown in Fig. 15 A ratio of 4 appears optimal, yielding a 
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field-enhancement factor of ~ 10. 



Figure 13. Planes showing a cross section at which potential (left) 
and electric field (right) values are calculated. One can see charge 
accumulation around the tips. 




Figure 14. Cross-section of the potential value - left, and the 
electric field strength- middle and right. Values are represented by 
color. The two-dimensional cross section is protruded in the direction 
perpendicular to the plane proportionally to the potential/field value 
at a given point. In this particular example, the spacing between the 
geometrical apexes is 4 times smaller than their amplitude. Note how 
quickly the edge effects disappear as one approaches the middle of the 
cross-section plane. 

What we conclude from the results in Fig. [15] is that optimal 
field enhancement is achieved at a case when the amplitude is 4 times 
larger than spacing between spikes. Starting from a perfectly flat 
surface, which has field enhancement of 1, the introduction of spikes 
(corrugation) starts yielding field enhancement. But as soon as the 
density is too high, field screening onsets and begins to counteract 
the field enhancement. At amplitude spacing ratio of 10, the system 
is almost equivalent to a smooth, flat electrode. These results show 
the strong, non-trivial dependence on the geometry of the system 
and how optimization needs to be carefully modeled. Since Maxwell's 
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Figure 15. Field enhancement values for different amplitude/spacing 
ratios for the egg-carton geometry. 

equations are scale- invariant, the geometry of the egg-carton used in 
this simulation applies equally at the nano scale, as well as the macro- 
scale without changing our conclusions. However, even minute changes 
in the details of the geometry (adding more spikes, varying height, 
etc.) will yield different behavior on optimal parameters for field 
enhancement. Therefore, the best option in optimizing such systems 
is an accurate model of the geometry. With a parallel version of the 
Robin Hood method, modeling of such geometries is now possible. 

4.6. Ion emitters 

Our second example is taken from the large field of research of ion 
emitters and gas detectors. In these systems, a neutral atom or 
a molecule is ionized once it approaches the electrode and is then 
accelerated in the electric field. Ion emitters have a wide variety of 
uses, such as compact neutron generators for oil well logging [25]. 
Even though average ionization energies are relatively small - few tens 
of eVs~, the ionization process is strongly determined by the shape 
of the potential barrier, i.e. by the electric field strength. Typical 
field strengths are on the order of 10^ — 10^'' V/m, and achieving 
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them with smooth, compact electrodes is difficult. Such difficulties 
can be overcome by having a geometry that can significantly enhance 
the electric field. Recently such devices have been producing using 
different nano structures, in particular using carbon nanotubes. Again 
for this study we have chosen a geometry that resembles experimental 
efforts, such as shown in Fig. 
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Figure 16. Geometry of the nano-tube array. On one half the actual 
mesh of triangles used in calculation is shown. 

Following closely the previous example, we calculate different 
arrangements of the nanotube array by varying both the radius of the 
nanotube and the inter-spacing between tubes. Once the surface charge 
densities are obtained, we calculate the potential and field values in a 



cross section plane, as shown in Fig. 17 



Having all the data for the field enhancement as was done for 



the egg-carton example, we obtain the results shown in Fig. 19 
The maximum field enhancement factor found using this nanotube 
structure can be as high as x64 (for a nanotube with radius-to-inter 
spacing value of 500) 



The results for field enhancement in Fig. 19 once again demon- 
strate that one isolated nanotube yields largest field enhancement. A 
higher packing density creates a screening effect sets which reduces this 
enhancement. In order to optimize the current per surface area, one 
should couple these results with a model of current strength depen- 
dence on the electric field strength. This approach is beyond the scope 
of this paper, but the interested reader should consult [30]. Again, in 
more realistic cases, the nanotubes may possess variations in orienta- 
tion and size that may alter the results shown here. None of these 
features poses a problem in principle, as long as one can provide a 3D 
mesh that describes the desired geometry with sufficient accuracy. 
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Figure 17. Cross section plane showing the electric field values. On 
the nano-tubes the charge density values are shown by means of color. 
Red color means large charge density, and one can see that the last row 
of nano-tubes is slightly more charged than the internal ones, having 
as the maximally charged nano-tube the one on the corner - this is also 
reflected in the electric field values in the cross-section plane. 

5. SUMMARY AND FUTURE WORK 

In conclusion, we have demonstrated the versatility of the Robin 
Hood technique across a variety of scales. In the case of large 
scale structures, such as the KATRIN detector, accuracy can be 
retained by introducing a high degree of segmentation and taking 
advantage of the low memory footprint used by the technique. In 
the case of micro-scale structures, we have shown that Robin Hood 
is capable of calculating large periodic arrays of complex structures. 
Microstructuring of electrodes is playing increasingly important role in 
many current applications ^34] • Simulations such as Robin Hood allow 
to move beyond simple approximations and provide more accurate 
optimizations. 

In principle, the boundary element method in conjunction with 
the Robin Hood technique can be made to solve any positive-definite 
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Figure 18. Electric field value in a cross-section plane shown from 
two different angles for the nanotube geometry. The cross-section 
plane is flat and it is being extruded in the perpendicular direction 
proportionally to the field value. The larger field enhancement at the 
very end of the array is due to finite size effects and are not taken into 
account in the the field enhancement calculation. 

matrix that exibits linear dependence. As such, we believe the 
technique can be utilized in other contexts, such as time-varying 
electromagnetic fields and linearized gravity. Such extensions are the 
subject of further inquiry. 



ACKNOWLEDGMENT 



J. A. Formaggio is supported by the United States Department of 
Energy under Grant No. DE-FG02-06ER-41420. The authors wish to 
thank the KATRIN collaboration for their support. 



Appendix A: Potential from Triangular Surfaces 

Let a triangular surface be defined by a corner Pq, the lengths of its 
sides a and b, and the unit vectors defining its sides ni and n2 (see 
Fig. [20| ). During the computation of the electric potential and field 

at a field point P, it is convenient to transform into a local coordinate 
frame §1 where the triangle lies in the x-y plane and the field point 
lies along the ^-clXlS, clS depicted in Figure ^Ol In the local coordinate 
system, the parameters necessary for the potential calculation are: 



Z = (P - Po) • ^3, 
as recommended in [35) . 



(47) 
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Figure 19. Field enhancement values depending on the radius and 
inter-spacing of the nano-tube array. The red X's show impossible 
configurations where nano-tubes would intersect one other. Notice 
that simple scaling is not possible in this case since the height of the 
nano-tube is kept constant in all cases. Units are arbitrary. 



X2 = xi + a, 

xy, = xi + (ni • n2) • 6, 



(48) 

(49) 
(50) 
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Figure 20. A triangular surface defined by a corner point Pq, the 
length of the longest side a and corresponding height b, and the unit 
vectors in the directions of the sides connected to Pq, labeled ni and 
n2 (ni always points in the direction of a). The field point is defined 
as P, with local coordinates (0,0,2;). The corners of the triangle are 
recast into local coordinates to facilitate integration. 



2/1 

2/2 



z 
\z\ 

2/1 



(p„-p) 

FT • 



n. 



2i 



(51) 
(52) 



where 713 = Hi x no- and is the unit normal to ni in the direction 
ofn2 (4 = ^^5^445^). 

^ ^ ^ |n2 — (n2-ni)ni| / 

The analytic calculation of the potential from this triangular 
surface (assuming a constant charge density a) can be found by directly 
evaluating the following integral: 



V = 



a 



where oi 



47reo 

0.2 



1 



• dx ■ dy, 



ai+biy Va;2 + y2 _^ ^2 

^^a^ and h 



(53) 



(y2-yi) ' ^ (yz-yi) ' ?/2-j/l ^ j/2-j/l 

We first manipulate the equation into a more soluble form: 



V 



o 

47reo jj/i 



/ / 



cosh (n) • 
sinh^ (ti) + 1 



where u\ = sinh 
identity 



-1 



and U2 = sinh ^ ^ 



cosh^(x) — sinh^(x 



1, 



(54) 

I . Using the 
(55) 
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Equation [54] becomes 

J U\ 

^ _o_ ry2 J 

47reo Jyi ^ 



du 



sinh 



(56) 



After dividing by z to make the integral dimensionless, Equation 
I becomes 



^ 47reo ^ 



au • smli ' 

Ul 

du ■ sinh 

Ul 



-1 f "'1 + ^1"" 



(57) 



where a\ = a'n = —, u = ^ and Ui = ^. 

If we define the following indefinite integral 



Ii{a,b, u) 



sinh 



-1 



a + 6u 



du, 



(58) 



Equation 57 can be rewritten as 
V 



4^ • z- [h{a2,b2,U2) - Ii{a2,b2,ui)- 
Ii{a[,bi,U2) + Ii{a[,bi,ui)] 



(59) 



5.1. Solution for Ji (a, 6, u) 

The solution to Ii{a,b,u) presented here is based on the technique 
outlined in [35]. Ii is first integrated by parts: 



h{a,b,u) = Fi{a,b,u) + h{a,b,u) - h{a,b,u), 



where 



Fi{a,b,u) = tt • sinh ^ 



a + bu 



73(0,6, u) 
14(0, b, u) 



au^ • du 

C/5(u2 + l) 

bu ■ du 

miu"^ + 1) 



(60) 

(61) 
(62) 
(63) 
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and the substitution 

[/= ((l + a2) + 2(a6)u + (l + 62)n2) (64) 

has been made merely to clarify the equations. By noting that Is can 
be rewritten as 

13(0, b,u) = a- I [ - [ — : 1 ' 

we can move the second term into 14, leaving us with 

/i(a, 6, u) = Fi{a, b, u) + ^{a, 6, u) — Ii{a, b, u), (66) 

where 

l3{a,b,u) = / (67) 
~ , , , [ ibu + a) ■ du , , 

By affecting a change in variables, 

a' = l + b\ 
y = 2ab, 

c = 1 + (69) 
I3 can be cast into a form readily found in tables of integrals: 
I3 f du 



- , , =. (70) 

a J ya'u^ + b'u + d 

The general solution to this integral for a' > is taken from ^36j to be 

/■ du 
J ^/a'u^ + Uu + d 

—= In (la'u + b' + 2^'^Ja!u^ + 6'n + d\ . (71) 

Recast into our variables of interest, we arrive at an equation for /s: 
l3{a,b,u) = • In (2 (6(a + 6u) + (72) 

Vi^TlVa'^ + "^abu + (6^ + 1) + l' 
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Integral I4 can be converted into a form with an analytic solution 
by performing the following change of variables: 



a 

b' 

c 

h 

k 



The recast integral is then 
h-- 



= l + 

= ab, 
= l+a^ 
= b, 
= a. 

{hu + k)du 



;2 + l)yja'v? + 2Vu + d 



(73) 



(74) 



and can be solved using a technique outlined in |37j by making the 
substitutions 

n -,du^--, — fdt (75) 

t — X 



in Eq. 74 and recovering the form 
h-- 



{t - xy 

{It + m)dt 



(t2 + 1)^0^2 + 2/3t + 7' 
This can be accomplished by making the following substitutions: 



(76) 



I = -hX - k 

m = kX — h 

a = X^a' + 2b' X + c' = 

= (l + 62)A2 + 2(a6)A + (l + a2), 
-b'X^ + [a' -c')X + b\ 
a' - 2b' X + c'A2 = 
= (l + 62)_2(a5)A + (l + a2)A2. 

Now, if we let A equal one of the roots of /?, 



—bX — a, 
aX — b, 



7 



A 



{a' - c') ± V(c' - a')'' + Ab''' 



or 



26' 

a 



(77) 

(78) 
(79) 



we can remove /3 from Equation 76 and split the integral into two parts: 



ltdt 



+ 



(t2 + l)y^^^ J (f2 + l)y^^2^ 



mdt 



(80) 
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The second integral can be put in the same form as the first if we affect 
another change of variables, t — t- and take t > 0, s > 0: 



ltdt 



msds 



+ a 

Solutions to integrals of this form are defined as 



(81) 



tdt 



tan' 



Using Equation 82 and taking the more negative root of /3 (A 
we arrive at a closed form solution for I/^: 

.2 1 



l4{a,b,u) = +b) ■ 

b V7 - a 



tan 




(82) 



(83) 



where 



a 

7 
t 



1 + 



(a2 + ^2) („2 + 52 + 1) 



62 



bu + a 



au 



(84) 



For computation, it is often beneficial to use an analytic form for 
14(0,6,^2) — 74(0, 6, Ml) to avoid intermediate computations that may 
be complex: 



h{a,b,U2) - h{a,b,ui) 



(85) 



+ 6 



1 



V7 -a, 



tan 



V7 - " • (^y^7*2 + «--\/7*i + a 
^ (7 - a) + y^Tti + a • ^Tti + « y 



It is also worth noting that, if the integration crosses a divergence in 
the integrand of (if ui < ^ < U2 or ui > ^ > U2), it is necessary to 
perform a branch cut on the interval of integration. 
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When z = 0, the above formahsm cannot be appHed, since it 
would involve a division by zero. Instead, we start with a modified 
version of Equation [56{ 



Ul 



^ dv 



sinh 



1/02 + b2y 



\y\ 



sinh 



(ai + biy 



V \y\ 

These two integrals can be made more soluble by integrating by parts: 

-1 f 0-2 + b2y" 



V 



4-Keo 



y sinh 



\y\ 



sinh 



I ( a i + biy 

\y\ 



y2 



+ 



y2 
yi 

y2 
yi 



dy 



(1 + bl)y^ + {2a2b2)y + 4 



dy 



(87) 

'il + bl)y^ + i2aibi)y + al 
The integrals in Equation l87] can be performed using the identity in 



Equation 71 
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